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■ We investigate the effects of the regularization procedure used in the J-Matrix method. We show that 

it influences the convergence, and propose an alternative regularization approach. We explicitly perform 
some model calculations to demonstrate the improvement. 

1 Introduction 

The J-Matrix (JM) method has become a popular tool in solving quantum mechanical scattering problems 
in atomic, nuclear and molecular systems. It is based on a square-integrable form of the wave function. 
The basis used for the representation of the wave function has to reduce the reference Hamiltonian, which 
is responsible for the asymptotic behavior of the solutions, to a tridiagonal or Jacobi matrix form. The 
scattering boundary conditions can then be expressed in terms of the expansion coefficients of the wave 
function. The Schrodinger equation has to be solved by matching the interior to the asymptotic region, 
which can now be done on the expansion coefficients. This leads to a set of linear equations, the solutions 
of which are the expansion coefficients of the interior representation, and the scattering parameters (phase 
shift or S-matrix). 

The JM approach proved to be succesful for 2-, 3- and many-particle systems in atomic and molecular 
DEI 02 il E|, and nuclear physics 00013111311111111131 

Despite the successes, a slow convergence as a function of the number of basis functions in the internal 
region leads to the need of large matrices. The calculation of matrix elements of the potential energy operator 
usually represents the main computational cost in a realistic many-body application. An improvement on 
the method leading to faster convergence with a smaller size of the Hamiltonian matrix is therefore necessary. 
There have already been attempts, based on asymptotic properties of the potential |14| . that led to sometimes 
dramatic reductions of the size of the modelspace. Other possibilities lie in computational methods for faster 
evaluation of potential matrix elements. Still a number of problems remain to solve the convergence issue. 

In this paper we analyze the behavior of the irregular (or Neumann-like) solution of the free-motion 
Schrodinger equation, and the corresponding expansion coefficients. Because of the square-integrable rep- 
resentation, a regularizing boundary condition has to be introduced to solve this equation. The standard 
regularization considered in the JM approach |15l 1161 II Tj will be shown to put a limit on the minimal size 
of the modelspace, and thus to have an important impact on the convergence. An analysis of the regu- 
larized asymptotic solution leads to a new regularization procedure that improves the size limitation. The 
enhancement in convergence of this procedure will be explicitly demonstrated on a number of potential 
problems. 
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We consider a JM formulation based on the oscillator basis, and a free-motion Hamiltonian. The analysis 
and the new regularization procedure can be easily extended to other bases, with analogous conclusions. 



2 J-Matrix Methods 

A Schrodinger equation for a continuum wave function with a spherically symmetric, non-coulombic potential 

h 2 ( 1 d ( 2 d \ 1(1 + 1) 



2m \r 2 dr\ dr ) h 2 r 2 ^ ^ J ^ ^ ^ ^ 

must have a solution that is matched asymptotically with the free-space Bessel and Neumann functions 

*;(r^oo)^ y^j ; (fcr)-tan<5;(fc) J^-m(kr) (2) 

where k 2 = ^SL^.This match of a solution in the interaction region, where the effect of the potential is felt, 
with the asymptotic reference states determines the phase shift at momentum k corresponding to energy 
E. We use the traditional spherical Bessel and spherical Neumann function definitions |181 IT§] with the 
delta-function normalization. 

The JM method reformulates this problem into an algebraic setting by the use of a basis functions. In 
our case these are the harmonic oscillator functions 

*««(»■ I h ) = (l) l L l n +1/2 ((^) 2 ) ex P (-|(J) a ) (3) 



N nl 



I 2T(n + 1) 
Y{n + l + i/2) 



where the L"(r) are Laguerre polynomials, N n i is a norm- factor and b is the oscillator length. We take a 
functional notation in which the variable of the function is separated from the parameter by a "|"; it will be 
used througout the paper. While it is straightforward to turn equation using the superposition ansatz 

oo 

*i{r\k) = X) Cni (*, (r 1 6) (4) 

into a matrix equation, more care must be taken when mapping the boundary condition J2J). The Bessel and 
Neumann functions are the solutions of the free-space reference Hamiltonian. 
The regular Bessel solution 

(T}-£0*f(r|fe)=O (5) 
*f(r|*)=y|M*r) (6) 
presents no problems and can be reconstructed as the solution of 

oo 

X ^nl(bm-E\^ ml (b))C^(k,b) = (7) 



rn—0 
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where the parameter of the basis states in brackets is indicated as a variable notation; this notation is used 
throughout the paper. This is a three-term recurrence relation that can be solved explicitly [1511161 IT?] : 

CZ (M) = N„,b"> (tt)'exp (-jit*) 2 ) L'* 1 ' 2 ((M>) 2 ) 

4fP|i ( ' ,, '*T !l ' F ' (-"' , + lot) (8) 

whre the the \F\ stands for Kummer's function [201 - F° r l ar g e n , the asymptotic behavior of the coefficients 
is given by 

C% (k, b) -> b^2R~^f (R nl | fc) (9) 



where the i? n ; = by/An + 21 + 3 are the oscillator turning points. 
The Neumann solution is irregular 

(T, - £!)*f(r|fc)=0 (10) 

*r(r|*)=yfnj(*r) (U) 

and cannot be mapped directly because the oscillator states cannot represent the singular behavior at the 
origin. This is not a fundamental problem because we only require these reference solutions to formulate the 
asymptotic conditions i.e. for r — ► oo. An equivalent solution, i.e. one that behaves asymptotically as the 
Neumann function but is regular at the origin, will do just as well. In some applications one uses a simple 
cutoff procedure. In the JM method it is of course important that the procedure translates easily to the 
algebraic representation of the problem. One takes the solution of 



(V-E)*? (r\k,b) = P$ ol (r\b) (12) 
tyf (r | k, b) — > \J~—ni{kr) for r — > oo (13) 

The regularized Neumann function can be determined from this equation using the Green's function technique 
[15|.|17| (see also next section). The coefficient (3 is fixed by the boundary condition. This regularized 
Neumann function can be reconstructed as a solution of the equation 

oo 

]T (<S>m(b)\T l ~E\<P ml (b))C? nl (k,b)=f3 S n , (14) 

m=0 

It is a three-term recurrence relation with a simple inhomogeneous contribution and can again be solved 
explicitly 



CZ (fc, 6) = y (fcb)-'- 1 exp (-^) iFx[-n-l--,-l + -;{kbYj (15) 

For large n, the asymptotic behavior is 

Cni(k, b) -> b^2R^f (r | A, 6) - by/tR^V? (R nl \ k) (16) 

We use these asymptotic values {JJJ and iflfi)) to compute the C^{k,b) and C^{k,b) coefficients by seeding 
the recurrence relations at some large index value, and recurring back towards small n. This turns out to 
be faster and numerically more efficient than directly evaluating the analytical expressions. 
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Thus we are led to a formulation of the scattering problem where we represent the solution as follows. 
The wave function has two essential components. In the internal region the effects of the potential are felt. 
We assume it has a finite range and that this component can be approximated by a finite combination of the 
L 2 basis states. In the asymptotic region the solution must match the boundary condition which we now 
express with the Bessel and regularized Neumann functions. 

z (r\k) = #f (r|fc) + *f (r\k) - tan^f (r\k, b) (17) 
*f(r|fc) for r -► oo (18) 

In the algebraic representation this becomes 

C nl {k, b) = C&(k, b) + C*(k, b) - tan (5, C%(k, b) (19) 
C^(k,b) ->0 for ra-> oo (20) 

The equation for the internal wavefunction and the phase shift is obtained by inserting (|17f) into the 
Schrodinger equation, projected onto the basis set: 

^$ Ili (6)|T / + ^-S|*f(/c) + *f(A:)-tan(5 ; *'f r (fc,6)^ = for all n (21) 
Taking advantage of the properties of the asymptotic reference states and (|12|l this simplifies to: 

(#„l(&)|T, +V-E\9l(k)) - tan<5, (f36 n0 + ($m(b)\V\*? (k, &))) 

= -<*„,(&)M*f(*)) (22) 

In the approach of [151 117| the algebraic Schrodinger equation is approximated by assuming that the in- 
teraction matrix for potential V can be truncated at some large n — N. This N defines a sharp bound- 
ary between the internal and asymptotic regions in the coefficient space. The resulting problem can be 
solved with C^ ; (fc,6) = for n > N . One has an N + 1 by N + 1 matrix equation in the unknowns 
{Pqv^i,V ■ ■ ■ i^n — i; ' tan Si}: 

AT-l 

£ (*ni(b)|T, + V - E\<f> ml (b)) CUk,b) 

m=0 

-Hm6i(k) (pS n0 + ^ni(b)\V\^ ml (b))C^(k,b)) 

\ m=Q / 

N—l 

= -J2(^i(b)\V\^ ml (b))C^(k,b) (23) 
Convergence of the results is achieved by extending the interaction region i.e. increasing N. 

3 The Regularized Neumann function 

The regularized Neumann function is a solution of the inhomogeneous differential equation (see details for 
the definition of *f (r\k, b) in HSHHini): 

(Tt-E)*? (r\k,b)=(3$ ol (r\b) (24) 
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where 

H 2 1 

^""^fc^o, (jfe|l/6)' ^ 

N 

The integral representation of ^ L (r\k, b) in coordinate space can be written as: 

oo 

*f (r I k,b) = 13 J G(r, r') $ /(r' | b)dr' (26) 
o 

where G(r, r') is the Green's function, explicitized in this case as: 

G{r y ) = d^ jl{kr<)ni{kr> y (27 ) 

In what follows we will restrict ourselves to the case of zero angular momentum. The regularized Neu- 
mann function can then be calculated analytically starting from the expression above, by straightforward 
integration of a combination of sines, cosines and gaussians. One finds 



< (r\k,b) 



= — [ cxp(— ikr) erf 
2-irkr 

+ exp(ifcr) erf 



ikb 

71 



/ kr 
\V%kb 



kr ikb 
V2kb 



(28) 



From this expression, we can see that the oscillator length has a strong impact on the behavior of (r | k, b). 
Indeed, when kb tends to zero, and if we take only the leading term in the series approximation of the complex 
error function for small complex argument, then 



< (r|M) 



2?rkr 



cos(kr) ■ erf 



/ kr \ 
\V2kb) 



(29) 



£ (r | k) • erf 



V2b 



N , 



Thus the larger is kr/kb — r/b, the closer the regularized function \& (r\k, b) is to the original Neumann 



function ^fff (f\k). For instance, when r > 56, the functions ^ (r\k, b) and (r|fc) are equal within single 
precision computation. When kb becomes large then 



—N 



2 1 



o(rM)«-\/-^cos(fcr) + --^ 



kb 



= ¥ N) (r\k) + - 



kb 



7T (kr) 2 + (kb)' 



HkW exp 

f l(kr] 
■ cxp < — - 



2 (kb) 



2(kb) 2 2 m 



2 --Akbf 



(30) 



In this case kr should be large compared to kb to suppress the second term in this expression, so that the 
function ^ (r\k,b) coincides with ^>q (f\k). The larger is kb, the larger the coordinate r has to be to reduce 

the difference between \P (r\k, b) and *S?q ( r \k)- For values of r < b , the regularized Neumann function can 
be represented as 



tf(,|M)— ^ 



— 2exp 



■2skberi(^ 
VV2 



+ 



(31) 
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by taking the leading term of the Taylor expansion of (|28[1 . Unlike the appearance, this expression is real 
because the error function is imaginary at an imaginary argument. One remarks (confirmed later by figure 
^) that in the limit for kb — ► one has 

—N 1 

* (0\k,b) ~ -— — ► -oo forfcb^O (32) 

/CO 

whereas when kb — > oo, the function behaves as: 

<(0|fc, fc) * 6Xp(( ff /2) _> +oo for kb - oo (33) 
fco 

The second situation in particular may be the source of numerical difficulties in actual calculations. Perhaps 
more to the point for the JM method is the fact that this behavior is reflected in the corresponding expansion 
coefficients C„ (k,b) for kb <C 1 and for kb 3> 1. Those coefficients then become very large at small values 
of n. From the analytical expressions ljl5|l for C Q N (k,b), and using the expressions in |20j for the limiting 
behavior of Kummer's function, one finds for instance: 



C ^(fc, b) = O ( ^ ) — ► +oo for kb -> (34) 



and 



exp ( ~ (Wc) 2 



C ^(fc, 6)=0 | ^ '- | — ► +oo for kb -» oo (35) 



Care is needed in these situations to avoid numerical difficulties associated with the large values of these 
coefficients.. 

In figure ^ w e display the behavior of the ^f (r\k,b) function for k = 1 and for values of the oscillator 
radii b = 0.1 up to b = 4.0 fm. One clearly notices that for small values of b the functions (r\k,b) and 
(r\k) are very close to each other in the whole range of r, with the exception of small r. The larger is 

b (or better, in general, kb), the larger is the value of the coordinate r where the functions ^ ( r l^i an d 
(r\k) are equal within required numerical precision. 

The present analysis has shown that the combined parameters for the momentum fc, and the oscillator 
length b (more specifically the product kb) have a great impact on the behavior of the function '>S l (r\k,b). 
For increasing b, given a fixed k, there is an increasing region of the coordinate r where the regularized and 
original Neumann functions differ considerably. Small values of b thus seem preferable for obtaining a faster 
convergence in the JM method. 



4 The J-Matrix method revisited 

In this section we will define new expansion coefficients for the regularized Neumann functions in order to 
increase the convergence rate of the JM method. For this aim we rewrite the equation (|24() in the following 

way 

(T, - E) V? (r\k, b ) = O $ o; (r\b ) (36) 

We have chosen an oscillator state $oz ( r \bo) for the regularization procedure that has a different oscillator 
radius bo than the oscillator radius b of the basis for the expansion {3J, and a corresponding /3q. The 
expansion coefficients of the (new) regularized Neumann function (|36(l 

C»> b (k,b ) = {$ni(b)\*?(k,b )) (37) 



G 



Figure 1: Neumann and regularized functions, for / = 0, k = 1 and various oscillator radii b. (Wj 
for b = 0.1 fm coincides with on this scale) 
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now satisfy the following set of algebraic inhomogeneous equations 



J2 ($ni(b)\T l -E\$ m i(b))Cy(k,b )=f3 {<Z> nl (b)\$ ol (b )) (38) 

m— 

i.e. a three-term inhomogeneous recurrence relation, with a source containing the overlap \*&oi{bo)). 
The latter can be easily calculated (see Appendix 0): 



<-<«i-<w-ii£,)7^(fc3)" 



where j = 1 + 3/2. This overlap decreases, for b > bo, for large values of n as 




(<m&) i $«(6o>« i^j Jf9y ex p(-( 2 "+^)|) ( 4 °) 

The larger is the ratio &o/fr, the faster the overlap tends to zero. When the ratio is small, the overlap slowly 
approaches to zero. In the region of n, where the overlap is negligibly small, the expansion coefficients 
C^ l ' b (k, bo) coincide with (k, b), defined by the formula fT5|l . 

The system of linear equations (|38J) can now be solved in the following way. We start the three-term 
recurrence relation in a region of n where the overlap is negligibly small, and take the original asymptotic 
^ni (fc>&) values as seeding values, then solve (|38f) for C^ h (k,bo) by stepping down the recurrence towards 
n = 0. 

Thus we are led to a formulation of the scattering problem where we represent the solution as follows. 
The wave function has two essential components. In the internal region the effects of the potential are felt. 
We assume it has a finite range and that this component can be approximated by a finite combination of the 
L 2 basis states. In the asymptotic region the solution must match the boundary condition which we now 
express with the Bessel and new regularized Neumann functions. 

tf, (r\k) = *f(r|fc) +*f (r\k) -tan$(fc)~#f (r\k,b ) (41) 
$j(r|fc) ^ for r -> oo (42) 

In the algebraic representation this becomes 

C nl {k,b)=C I nl {k,b) + C% l {k,b)-ta,n5 l {k) C% b (k,b ) (43) 
C^(k,b) ->0 for n -> oo (44) 

The equation for the internal wavefunction and the phase shift is obtained by inserting (|17Jl into the 
Schrodinger equation, projected onto the basis set: 

( $„K&)|r/ + y-S|*f(fc) + ^f(fc)-tan^Wf (fc,&o)) =0 (45) 
Taking advantage of the properties of the asymptotic reference states (JjJJ and ()12|l this simplifies to: 



(^(bm+v-Ei^ik)) 

i(b] 

= -(<M6)|V|*f(*0) (46) 



tantf, f/3 ($nz(&)|*f (k,b )) + U nl (b)\V\V? (Mo) 
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Figure 2: s-wave phase shift for the square- well potential, obtained with b = 3 fm and varying renormalization 
widths b 



In the approach of [151 I17| the algebraic Schrodingcr equation is approximated by assuming that the 
interaction matrix for potential V can be truncated at some large n = N. This N defines a sharp 
boundary between the internal and asymptotic regions in the coefficient space. The resulting problem 
can be solved with C^;(fc,6) = for n > N . One has an N + 1 by N + 1 matrix equation in the 
unknowns { C^, C[ C I N _ lll tan Si}: 

N-l 

( Ti + V-E\ $ m; (6)> C^ik, b) 
-tnnSt Uc^ 6 (fc,6o) + E ^ni(b)\V \$ ml (b)) C% b (k, b ) j 

\ m=0 / 

N-l 

= -J2 (®m(b)\V\$ ml (b)) C^(k,b) (47) 

m— 

Convergence of the results is achieved by extending the interaction region i.e. increasing N. 

5 Some examples 

In this section we present some detailed results for radial 1-dimensional model potentials. We compare the 
application of the traditional JM rcgularization procedure to the one introduced in section 0] for Gauss, 
exponential, Yukawa and square-well potentials. We choose comparable parameters, in particular a depth 
of Vo = —80 MeV and width a = 1.0 fm, for all of these potentials. 

All calculations are made in the standard JM approach, i.e. truncating the potential matrix beyond 
the boundary condition matching point N ^ 1151 121j . For the traditional regularization we consider the 
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Figure 3: s-wave phase shift for the Yukawa potential, obtained with b — 3 fm and varying renormalization 
widths bo 



standard C^i(k,b) and C^i(k,b) asymptotic expansion forms, whereas for the new regularization (J-Matrix 
Regularized (JMR) method) the C^(k,b) and C^' b (k,bo) asymptotic forms are used. We have set the 
matching point to N = 25 in all cases. To obtain a fair comparison of the convergence properties for the 
two regularization properties, we compare the s-wave phase shifts to each other and to the "exact" result, 
obtained with the Variable Phase approach |22l l2"3*| . 

In figures [3 and |21 we show how the solutions obtained with the new regularization (JMR) depend on 
its defining value bo. We have done this for both the Yukawa an square- well potentials, as it is on these 
potentials that the effects are more noticeable. It is also seen that the optimal value of bo strongly depends 
on the functional form of the potential, more particularly on the behavior at the origin. 

In figures 01 El and the comparison of the JM, JMR and VPA s-wave phase shifts for all four model 
potentials is displayed, for three different values of the oscillator length b — 1.0, 2.0 and 3.0 fm. For 
uniformity the value for the regularization parameter was chosen to be bo — 0.6 fm in all cases, as it is a near 
optimal value to limit the differences between the renormalized <J L (r\k,bo) an d true asymptotic Neumann 

( r \k) functions in coordinate representation (see figure^. 

We notice that, for all potentials, the results deviate more from the exact (VPA) results with increasing 
b. This is an indication that a larger value for the matching point N should be chosen for truly convergent 
results. The phase shifts obtained by JMR however remain much closer to the VPA results. The potential 
effects are more pronounced for the square- well (figureEJ) and Yukawa potential (figured). The regularization 
effect is in all cases seen to be more important for large values of b where the convergence problem becomes 
an important issue. 

It is also clear that convergence, both for JM and JMR, can be more easily achieved for the Gauss and 
exponential potentials than for the square-well and the Yukawa potentials. This undoubtedly has to do with 
the more regular behavior at the origin of the former potentials. The JMR results can be further improved 
if a more judicious choice is made for the regularization parameter bg. Indeed, from figures [3 and |31 it is 
seen that a value of bo around 0.4 fm (square-well) or 0.2 fm (Yukawa) would improve the results even more 
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Figure 4: Comparison of s-wave pahse shifts for JM and JMR for a Gauss potential (Vo = 80 MeV, a = 1.0 
fm). The matching point is TV = 25 in both calculations. The Regularization parameter is bo = 0.6 fm. 
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Figure 5: Comparison of s-wave pahse shifts for JM and JMR for an exponential potential (Vb = 80 MeV, 
a = 1.0 fm). The matching point is N — 25 in both calculations. The Regularization parameter is b = 0.6 
fm. 
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Figure 6: Comparison of s-wave pahse shifts for JM and JMR for a square-well potential (Vo = 80 MeV, 
a = 1.0 fm). The matching point is TV = 25 in both calculations. The Regularization parameter is b = 0.6 
fm. 
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Figure 7: Comparison of s-wave pahse shifts for JM and JMR for a Yukawa potential (Vb = 80 MeV, a = 1.0 
fm). The matching point is N = 25 in both calculations. The Regularization parameter is 6 = 0.6 fm. 
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drastically. 

It should also be clear that combining the new regularization method introduced here with other methods 
for improving convergence, e.g. using semiclassical potential considerations such as |14j . will carry the merits 
of both approaches, and can dramatically reduce the size of potential matrices. 



6 Conclusions 

We have investigated the regularization procedure of the JM method. We have shown that the differences 
between the regularized Neumann and Neumann functions can have a significant impact on the convergence 
of the JM method. We have proposed a new approach to the regularization. It is computationally only 
slightly more involved but it yields significant improvement in the convergence of the phase shifts obtained 
by the method. We have explicitly demonstrated the latter point through calculations on several model 
potentials. 



A Overlap 



To calculate the overlap between oscillator functions of different widths we use the technique of generating 
functions. It is known that the function 



<f>l(r\e,b) 



1 



+ 



-{1+3/2) /r\ J +3/2 



r(Z + 3/2) 6 3 / 2 

generates a complete set of the oscillator functions of width b 



exp 



1 1 — e /r\ 2 
2 1 + e U 



(48) 



-7fc 

$ i (r|e,&) = ^- r $„ i (r,b) 



(49) 



where 



Using the relation 



Til 



Mnl 



T(n+j) 



3 = 1 + 3/2 



=o = ®ni [r\ b) 



(50) 



(51) 



on matrix elements of the generator state, one can extract the oscillator function matrix elements for definite 
quantum number n. 

It is straightforward to calculate the overlap of two generating functions. One finds 



3,(e,&) | *, (?,&))= [V7£ 



where 



1 



(1+7)' 



(1-7) 
(1 + 7)' 



A = 1 — as + ae — ee 



(52) 



(53) 



Closer inspection of the expression will reveal that it is actually symmetric in the width parameters as it 
should be. By applying relation Q51JI to this overlap we obtain the required overlap of the oscillator functions: 



N n l Nnl 

n\ n\ 



IX (4X <<!,(.-./,) | <I.„7. U) 



de 



(54) 



J e=e=0 
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It can be represented in the following form: 

[$„,(!&) I *nl(| 6) 



3 „n+n T ( n +j) T ( n + j) 



n\nW(j)T(j) 



min(n,n) 

E 



n\n\T (J) 



k\ (n - k)\ (n - k)\T {k + j) 



I -a 2 



OL 



M i (^^" r(B+i)r( " + i) 



2 F 1 -n, -n;j;- 



n!n!r (j) V (j) 
(l-a»)\ 



For the particular case where one of the quantum numbers is zero one finds: 

, +a r(n + j) 



$nj(|&)|*nl(l 6)) =M,I («>/7) 



266 



n!r (j) 



T(n + j) /6 2 -6 2 



6 2 + 6 2 / V n!r 0") V b 2 + b 2 



(55) 



(56) 
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